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ABSTRACT 



A method for implementing cylindrical coordinates in the Athena magnet ohydrody- 
namics (MHD) code is described. The extension follows the approach of Athena^s orig- 
i nal developers an d has been designed to alter the existing Cartesian-coordinates code 
(IStone et al.l l2008l ) as minimally and transparently as possible. The numerical equa- 
tions in cylindrical coordinates are formulated to maintain consistency with constrained 
transport, a central feature of the Athena algorithm, while making use of previously im- 
plemented code modules such as the Riemann solvers. Angular- momentum transport, 
which is critical in astrophysical disk systems dominated by rotation, is treated carefully. 
We describe modifications for cylindrical coordinates of the higher-order spatial recon- 
struction and characteristic evolution steps as well as the finite- volume and constrained 
transport updates. Finally, we present a test suite of standard and novel problems in 
one-, two-, and three-dimensions designed to validate our algorithms and implementa- 
tion and to be of use to other code developers. The code is suitable for use in a wide 
variety of astrophysical applications and is freely available for download on the web. 



Subject headings: hydrodynamics - MHD - methods: numerical 



Introduction 



T he Athena code flGardiner fc Stond [20051, hereafter IGSOSI : iGardiner fc Stond [20081, hereafter 



GSOSi : [stone et al.[ [2008[) is a new, second-order Godunov code for solving the equations of ideal 
magnetohydrodynamics (MHD). Among its salient features are that it preserves the divergence- 
free constraint, V • B = 0, to within machine round-off error via unsplit evolution of the magnetic 
field, and that it employs fully conservative up dates of the MHD equa tions. This last feature 
distinguishes Athena from its predecessor, Zeus ( Stone fc Norman 1992al lbl). which also preserves 
the divergence-free constraint, but employs operator-split finite-difference methods. Athena has 
been extensively tested via both comparison to analytic solutions, and comparison to the results of 



askinner@astro.umd.edu 
^ostriker@astro. umd.edu 



- 2 - 



other numerical MHD codes. The code package is freely available to the community, and is highly 
portable and easily configurable, as it is self-contained and does not rely on outside libraries other 
than MPI for computation on multi-processor distributed memory platforms. 

The equations of ideal MHD consist of eight coupled partial differential equations, which are 
not analytically solvable in general, and fully three-dimensional numerical solutions can be quite 
costly. For many astrophysical systems of interest, however, the computational cost for certain 
problems can be reduced by exploiting geometric symmetry. For example, the high angular velocity 
of the plasma in accreting systems implies that most of the mass is confined within a disk. If the 
properties are statistically independent of azimuthal angle, 0, these disks can be studied using 
radial-vertical (R-z) models, and if vertical variations are of lesser importance, these disks can 
be studied using radial-azimuthal (R-c/)) models. The dynamical properties of winds and jets from 
astrophysical systems can also be analyzed using axisymmetric models. Exploiting symmetry in this 
way to reduce the effective dimension of the problem can greatly simplify the calculations involved 
and allow finer resolution when and where needed. In addition, for either reduced-dimensional or 
fully three-dimensional problems, using a curvilinear coordinate system for rotating, grid- aligned 
flow is superior for preservation of total angular momentum, and renders imposition of boundary 
conditions much simpler compared to the Cartesian-grid case. 

There are several other publicly available hirfj-resolution shock-c apturing!; codes for astrophys- 



i cal MHD in wide use, includ ing 1/^(7 (iTothI Il996ll. BA TS-R-US f|Powen et al.Ml999l). FLASH 



( Frvxell et al ] l200nl ). RAMSES ^^s^^oSh NIRVANA flzieglei][2004l l. and PLUTO dMignone et al 
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3), to name a few. Although these and other codes enjoy increasing popularity within the com- 



munity, as of this writing only VAC and PLUTO have the capability for MHD in curvilinear 
coordinates. 

In this paper, we describe our adaptation of Athena to support cylindrical geometry, and 
present a suite of tests designed to validate our algorithms and implementation. These tests include 
standard as well as novel problems, and may be of use to other code developers. A guiding principal 
of our approach is to alter the existing Athena code as minimally and as transparently as possible. 
This will involve a careful formulation of the MHD equations so that the finite-volume algorithm 
remains consistent with constrained transport, and so that the built-in Riemann solvers (as well as 
computation of wavespeeds and eigenfunctions) need not be changed. Finally, we pay particular 
attention to angular-momentum transport, which is critical in systems dominated by rotation. 

The plan of this paper is as follows: In ^ we describe the conservative system of mathematical 
equations that we shall solve, and in ^ we briefly outline the main steps used in Athena to evolve 
the system numerically. In we describe the projected primitive variable system used in the 
reconstruction step. In §^ and El we describe the modifications needed for cylindrical coordinates 
in the higher-order spatial reconstruction and characteristic evolution steps, respectively. In §^ 
andEl we describe the implementation in cylindrical coordinates of the finite volume and constrained 
transport updates, respectively, and then in ^ we summarize the steps of the whole algorithm in 
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detail. In ^TOt we present code verification tests and results, and we conclude in ^TTl 

Our version of the code, including the suite of test problems we have developed, is freely 
available for download on the Web. 



2. The Equations of MHD 



The coordinate-free conservative form of the equations of ideal MHD are: 

dtp + v-{pv) = 0, 

dt{pv) + V -{pw-BB + P*!) = -pV$, 
dtE + V ■[{E + P*)v-B{B-v)\ = -pv-V^, 
dtB + V -{vB-Bv) = 0. 



(la) 
(lb) 
(Ic) 
(Id) 



Here, p is the mass density, pv is the momentum density, B is the magnetic field vector, and I is 
the identity tensor. The total pressure is defined as P* = P + (S • S)/2, where P = nkT is the 
thermal pressure, and E is the total energy density defined as = e + p{v •v)/2-^{B- S)/2, where 
e is the internal energy density. An ideal gas equation of state P = (7 — l)e is assumed, where 7 
is the ratio of specific heats. As written, the equations have been scaled in such a way that the 
magnetic permeability is /i = 1 (for cgs units, B is replaced by B/y/^). Optionally, we can include 
a static gravitational potential $ = in equations (ITb]l and (ITc]l : the energy equation (ITc|l can 

also be generalized by including radiative heating and cooling terms. 

Ignoring terms on the right-hand sides, equations ^ can be summarized by the single evolution 
equation in "conservative" form: 

a^Q + V • F = 0, (2) 
where Q = Q{x,t) is the set of conserved quantities 

/ P \ 

pv 
E 

\ B ) 



Q = 



(3) 



and 



F = 



pv 
M 

{E + P*)v - B{v ■ B) 
J 



(4) 



is a structure whose components represent the (nonlinear) fluxes associated with the various com- 
ponents of Q. For added simplicity, we have defined the momentum and induction tensors: 



M. = pvv- BB + P*I, 



(5) 
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and 

J = vB - Bv. (6) 

Note that additional terms such as gravitational forces appearing on the right-hand sides of equa- 
tions ^ are treated separately as source terms, hence are not part of the conservative system in 
equation ([2]). 

In Cartesian coordinates, equation ([2|) can be expanded in a straightforward manner: 

dtQ + d^F, + dyFy + d,F, = 0, (7) 

where Fx, Fy, and Fz are one-dimensional vectors representing the fluxes of various components 
of Q in each orthogonal coordinate direction. 

In order to extend the existing algorithm in Athena to curvilinear coordinates, we introduce 
geometric scale factors and source terms arising from the covariant derivatives in the curved metric. 
Following this approach, the existing code can be made to support cylindrical geometry with only 
moderate adjustments, which is what we describe in this paper. 

We expand equation ([2]) according; to the form of the divergence operator in cylindrical coor- 



dinates (see e.g. iMihalas fc Mihalaslll984l ) when acting on a vector, 

V • = ^dn {Rvr) + ^dcfyV^ + d^Vz, (8) 

and when acting on a tensor, 

(V-T)^ = ^dR{RTRR) + ^d^T^R + dzTzR-^T^^ (9a) 

(V-T)^ = ^dR{RTR^) + ^d^T^^ + dzTz^ + ^T^R (9b) 

(V-T), = ^dR{RTRz) + ^d^T^z + dzTzz. (9c) 



The extra non-derivative terms in equations (I9al) and (I9bl) are the so-called geometric source terms, 
and represent "fictitious" forces, e.g. the centrifugal and Coriolis forces. Once source terms are 
introduced, the finite-volume updates are no longer fully conservative. As we shall next show, 
however, all but one of the geometric source terms can be eliminated from the equations. This 
remaining geometric source term, appearing in the radial momentum equation, is often balanced 
by gravity in realistic astrophysical problems. 

2.1. Continuity Equation 

Expanding the derivative operators in cylindrical coordinates in equation ([Tall , we have the 
continuity equation in conservative variable form: 

dtP+^dR{RpvR) + ^d^ipv^) + dzipvz) = 0. (10) 
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2.2. Momentum Equation 



For the momentum equation (jlbp . we have, in terms of the symmetric momentum tensor M 
(eq. El): 

dt{pvR) + ^dR{RMRR) + ^d^M^R + d,M,R = ^M^^ (11a) 
dt{pv4,) + ^dR{RMR^) + ^d^M^4, + d,M,^ = -^M^r (lib) 
dt{pv,) + ^dR{RMR,) + ^d4,M^, + d,M,, = 0. (11c) 



Mignone et alJ (|2007l ) note that the symmetric character of M ahows further simphfication of 



equation (jllbj) in cyhndrical coordinates, since 

(V • M)^ = -^dR{R^ Mr^) + ^d^M^^ + d,M,^. (12) 

This leads to the so-cahed "angular momentum-conserving form" of the (/)-momentum equation: 

dt{pRv^) + ^dR{R^MR^) + ^d^{RM^^) + d,{RM,^) = (13) 

Note that in this form, the conserved quantity is angular momentum and there is no source term. 
However, since the original Cartesian version of Athena makes use of linear momenta in the flux 
calculations and since we do not wish to alter those calculations, we rewrite this equation once 
more to obtain the ^-momentum equation that we shall use: 

dt{pv^) + ^dR{R^ Mr^) + ^d^M^^ + d,M,^ = 0. (14) 

This leaves the term M^j^^f^/R = {pv'^ — S| + P'^)/R appearing in the i?-momentum equation flllap 
as the only geometric source term in the cylindrical coordinate expansion of equation (llbl) . In 
practice, this source term is often (partially) balanced by the radial component of the gravitational 
source term 

-pV$ (15) 

for many astrophysical applications. 

2.3. Energy Equation 

For the total energy equation (fTcl) in conservative form, we have 

dtE + \dR[R{{E + P*)vR-BR{B-v))\ + ^d^{{E + P*)v^-B^{B-v)) 

+ d,{{E + P*)v,-B,{B-v)) = Q. (16) 

The gravitational source term for the energy equation is 

-pv-V^. (17) 
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2.4. Induction Equation 

Finally, for the induction equation (jldp . we have in terms of the antisymmetric induction 
tensor, J (eq. El): 

dtBR + ^d^J^R + d,J,R = (18a) 
dtB^ + ^dR{RjR^) + d,J,^ = ^ (18b) 
dtB, + ^dR{RjR,) + ^d4,J4„ = 0. (18c) 

It is important for the preservation of the divergence constraint that we avoid source terms in the 
magnetic fluxes. Thus, we rewrite the ^-induction equation (|18bp as 

dt[^)+^dR[R'-^)+d.[^)^^. (19) 

Note that in this form, the conserved quantity is B^j^/R and there is no source term. Equation (fT9]) 
has also modified the fluxes of B^j^, but since Athena does not actually use flux differences to evolve 
the magnetic fields (see this is not a serious problem. We will use a reduced form equivalent 
to equation (fTOll for B^jy in the reconstruction step (see 

dtB^ + OrJr^ + d,J,^ = 0. (20) 

In this form, we use the fluxes Jji^f) and J^^ originally appearing in equation (Il8bl) . not the modified 
fluxes JR(f)/R and Jzcfy/R appearing in equation (IT9]) . 



3. Overview of the Numerical Algorithm 

The algorithm used in Athena to evolve the system in equation ([2|) us es a Goduno v-type finite 
volume (FV) scheme. A simplified version of the algorithm presented in IStone et al. dioos.) IS as 
follows: 

1. Using cell-centered volume averages at time t"^, compute left and right (L/R) interface states 
with the reconstruct-evolve-average (REA) method based on the linearized one-dimensional 
evolution equations. 

2. Add the parallel components of source terms to the L/R states. 

3. Compute the first-order interface fluxes from the L/R states using an exact or approximate 
Riemann solver. 



4. Upda te the L/R states' magnetic fields using constrained transport (CT) (lEvans fc Hawley 



19881). 



5. Correct the L/R states' remaining non-magnetic variables with transverse flux gradients and 
the transverse components of source terms. 
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6. Compute the second-order interface fluxes from the corrected L/R states using the Riemann 
solver. 

7. Using the second-order fluxes, advance the interface mag netic fields to time t^+i = t + At 
with CT. 

8. Using the second-order fluxes, advance the remaining cell-centered quantities to time f^^^ 
with the FV method. 

9. Add the time- and volume-averaged source terms to the cell-centered quantities. 

10. Average the interface magnetic field components to obtain the cell-centered field components 
at time 

11. Compute a new timestep At based on the CFL condition and repeat steps ffB- ffTTI) until 

> tf. 



Currently, Athena includes a wide variety of non-linear Riemann solvers fs eelStone et al.ll2008l . 



for a complete list). In our tests, we use the solvers based on the HLL flux (IHarten et al.lll983l ) 



as well as Ro e's linearized method (iRodllQSlI ). Among the solvers based on the HLL flux ar e the 



HLL E solver (lEinfeldt et al.lll99ll ). which uses a single intermediate state, the HLLC solver (IToro 



19991). which extends the original HLLE solver by including a contact wave, and the HLLD solver 



^jniyoshi fc KusanQll2005l ). which extends the original HLLE solver by including both contact and 
Alfven waves. The HLLE solver has the advantage that it is simple and therefore faster than more 
accurate solvers such as Roe's method, and like all solvers based on the HLL fluxes, it is positive- 
definite for ID problems. However, since it neglects the contact wave, and additionally the Alfven 
waves for MHD, it is overly-diffusive for these waves. On the other hand, for hydrodynamics, the 
HLLC solver produces fluxes that are as accurate, if not better, than those produced by Roe's 
method, but at a considerably lower computational cost. For MHD, it has been shown that the 
HLLD solver is of comparable accuracy to the MHD extens ion of Roe's method for several tests 



using Athena^ although it is much faster (IStone et al.ll2008l ). The advantage of Roe's linearized 



method is that it includes all waves in a given problem, yielding less diffusive and, hence, more 
accurate results for intermediate waves that are neglected by the methods based on HLL fluxes, 
although for some values of the left and right states. Roe's method will fail to return positive density 
and/or pressure in the intermediate state(s). Finally, we reiterate that with the approach we have 
adopted, it is not necessary to make any changes to the computation of wavespeeds, eigenfunctions, 
or fluxes in any of these methods. 

Although no changes are required for the solution of the Riemann problem at interfaces, several 
changes are required in other parts of the Athena algorithm in order to accommodate non-Cartesian 
coordinates. In the next sections, we discuss the geometry-specific details of computing the L/R 
states (steps[l]|2l see §SlEj), the FV method (steps[8l|9l see and the incorporation of CT into the 
corner transport upwind (CTU) method of (Colellal fll99G ) (see Finally, we will recapitulate the 
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steps of the algorithm in greater detail and explain the computation of the new timestep (step [TTJ 
see 



4. The Linearized Evolution Equations 



In Athena, the left and right (L/R) interface states (the inputs to the Riemann solver) are 
computed using a modified form of the system in equation ([2j). The equations, written in primitive 
variable form, are projected in a single coordinate direction, and the resulting system is linearized 
and then evolved. The projection in the ^ -direc tion yields a system that can be obtained from 
the corresponding Cartesian projection (see IGS05 , §3.1) by making the substitution dy R~^d(p. 
However, the projection in the i?-direction differs more significantly as a result of geometric scale 
factors. 

For the projection in the i?-direction, we begin with the primitive variable form of equation ([2]), 
take d(f) = and dz = 0, expand the remaining i?-partials, and move the non-derivative terms to 
the right-hand side to obtain the system: 



dfW + AOrw = s. 



(21) 



where 

P 

vr 

W ^ Vz 
P 

B4 

is the vector of primitive variables, omitting the parallel component of the magnetic field. 



(22) 



Vr 


p 




















Vr 








1/p 


B4>/p 


B,/p 








Vr 








-Br/p 














Vr 








-Br/p 





IP 








vr 













-Br 








VR 













-Br 








VR 



(23) 



is the wave matrix, and s — smhd + <Sgrav + <Sgeom is the source term vector, a combination of the 
MHD source terms arising from the V • B constraint, gravity source terms from a static potential, 
and the geometric source terms inherent in the cylindrical coordinate system. As in the Cartesian 
version of Athena, the form of the MHD source terms differs slightly in the 2D and 3D cases (see 
] below), but the forms of the gravity and geometric source terms are independent of dimension. 



-9- 



The hyperbolic wave matrix, A, given in equation (l23l) . is linearized by taking it to be a 
constant function of the primitive variable state w at time f^. However, it is only indirectly 
accessed through the system of eigenvectors and eigenvalues of A (see ^ below). We write the 
projected equations in cylindrical coordinates using this specific form in order to make use of the 
eigensystem solution previously implemented in Athena. 

In the remainder of this section, we derive the cylindrical coordinate form of the primitive 
variable system given in equation (|2T]) . and in the process obtain the geometric source terms. 

4.1. Continuity Equation 

Expanding the derivative operators in cylindrical coordinates in equation (fTal) and projecting 
in the i?-direction, we have for the continuity equation in primitive variable form: 

dtp + pOrVr + vrOrp = -^pvR- (24) 

The left-hand side of equation ([24|) contains all the terms from equation ([2T1) . and the term on the 
right-hand side is the first component of the geometric source term vector, Sgeom- Furthermore, if 
we make the substitution R ^ x and ignore the source term, we recover the x-projection of the 
continuity equation in Cartesian coordinates. 



4.2. Momentum Equation 

For the momentum equation, we begin with the conservative form of equation (fTb|) and use 
the continuity equation and divergence-free constraint to eliminate terms and obtain: 

pdtv + p{v • V)v - {B • V)B + VP* = 0. (25) 

By explicitly enforcing V • S = here, we ensure that any numerical error in the divergence of the 
magnetic field can not influence the evolution of momentum during the reconstruction step. 

Next, we divide through by p, substitute P* = P + P^/2, project in the P-direction, expand 
the partials, and move the source terms to the right-hand side to obtain: 

dtVR^VRdRVR + \dRP^\B^dRB^^)^B,dRB, = ^{vl-\Bl), (26a) 
dtVct> + vrOrv^ - -^BrDrB^ = -^{vct^VR - ^B^Br), (26b) 
dtVz + vrOrVz - ^BrOrB^ = 0. (26c) 

Recall that the ^-momentum equation (|T3l) can be expressed in angular-momentum conserving 
form and thus avoid a geometric source term. However, we must include the source term on the 
right-hand side of equation (I26bl) in primitive variable form in order to preserve the specific structure 
of the coefficient matrix. A, on the left-hand side of equation (j21l) . Finally, the gravity source terms 
in the momentum equation are given by the components of — V$ in cylindrical coordinates. 
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4.3. Energy Equation 

We begin with the internal energy equation in coordinate-free form: 

dtP + v VP + -fPV v^O. (27) 

Then, projecting the equations in the i?-direction, expanding the partials, and moving the source 
term to the right-hand side, we obtain: 

dtP + vrOrP + -fPdRVR = -^jPvR- (28) 
In primitive form, there is no gravity source term in the energy equation. 



4.4. Induction Equation 



For the induction equation (fldl) , also written as dtB — V x (i; x S) = 0, we begin with the 
components given in equations fjl8al) . fjl8cl) . and f[20D . Moving terms proportional to R~^dR{RBji)^ 
R~^d(f)B(i)^ and dzB^ to the right-hand side, we obtain: 

dtBR + ^d^iv^Bn) - B^^d^VR 

+d,{v,BR) - B,d,VR = VR[^d^B^ + d,B,], (29a) 
dfB^ + dR{vRB^) - BrOrv^ 

+d,{v,B^) - B,d,v^ = v4^dR{RBR) + d,B,]-^v^BR, (29b) 
dtB, + ^dR[R{vRB,)] - BrOrv, 

+ ^d^{v^B,) - B^^d^v, = v,[^dR{RBR) + ^d^B^]. (29c) 

In 2D, the divergence- free constraint in cylindrical coordinates with = implies that the right- 
hand side of equation (|29cl) is identically zero. Applying the divergence constraint in equation (|29cl) . 
projecting in the i?-direction, and expanding the i?-partials on the left-hand side, we obtain: 



dfB^ + B^Orvr - 
dtB, + B^Orvr 



dtBR 

BrOrv^ + VrOrB^ 
■ BrOrVz + VrOrB, 



0, 

Vcf>^dR{RBR) 
-^vrB,. 



(30a) 
(30b) 
(30c) 



Note that the left-hand sides of equations f)30bD and f|30cl) exactly match the Cartesian form (jGSOSL 
eq. 30) if we make the substitution R x. Also, the divergence term on the right-hand side of 
equation (|30bl) matches the divergence term on the right-hand side of the Cartesian form, except 
that it appears in cylindrical coordinate form. However, the additional source terms —Vcj^Br/R in 
equation (|30bp and —vrBz/R in equation (I30cp . which vanish as i? ^ oc, are curvature-related 
terms that are unique to the system in cylindrical coordinates. 



In 3D, the cancellation of the V • B terms on the right-hand side of equation (I29cp is no longer 
possible. Instead, iGSOSi introduce an algorithm that adds a limited amount of the MHD source 
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terms from the transverse directions to the source terms in each sphtting direction. This is done 
in such a way that the overah induction equation is not altered and so that the sum of the MHD 
source terms is minimized. These constraints take the form of minmod hmiter f unction s that reduce 
to the underlying 2D algorithm in the limit of 2D, grid-aligned problems (see IGSOSI , §3). The 3D 
algorithm in cylindrical coordinates yields the system 

dfBR + {^d^{v^BR - B^vr) - VRLR^{dzB^)} 

Hdz{v,BR-B,VR)-VRLR,{^d^B^)} = 0, (31a) 
OtB^ + {dR{vRB^ - Brv^) - v^L^R^d^B^)} 
HdzivzB^ - B,v^) - v^L^,[^dR{RBR)]} = 0, (31b) 
dtB, + {^dR[R{vRB, - Brv,)] - v,L,R{^d^B^)} 

H^d^{y<pBz-B^v,)-v,L,4^dR{RBR)]} = 0, (31c) 

where 

LRcfy{dzBz) = mmmod{- ^dcfyB^, d^B^), (32a) 
L^dzB,) = minmodi- ^dR{RBR), d,B,), (32b) 
Lzni^dcfyB^) = minmod(-^ai^(i?Si^), ^d^B^). (32c) 

Note that the limiter Lij is only applied to the equation for Bi projected in the j-direction. We 
require that 

LRzi^d^B^) = -LR^{dzB^), (33a) 
L^z[^dR{RBR)] = -L^R{dM, (33b) 
L,4^dR{RBR)] = -L,R(^d^B^), (33c) 

so that the limit ers cancel pairwise when summed over all projections. The limiters defined in 
equations f[32l) and f[33l) are the same as in the 3D Cartesian formulae f|GS08l . eqs. 11,15), with 
d^B^ ^ R-^dR{RBR) and dyBy ^ R'^d^B^. 

Projecting equations ([HT1) in the i?-direction, expanding the partials (except for those appearing 
in the limiter functions), moving all of the source terms to the right-hand side, and using the 
properties of the minmod function together with the V • S = constraint, we obtain: 

dtBR = 0, (34a) 
dtB^^ B^Orvr^vrOrB^- BrOrv^ = v^miimiod[^dR{RBR), -^d^t^B^] 

-^v^Br, (34b) 
dtB^ + B^Orvr + vrOrB^ - BrOrVz = Vz mmmod[^dR{RBR), -d^B^] 

-^vrB,. (34c) 

Note that for the 2D case with 5^ = 0, the V • B constraint implies that the arguments of the 
minmod function in equation (I34bp are equal and that the minmod function in equation (I34cp 



evaluates to zero, so that we recover the 2D system in equations (l30]) . The minmod terms on the 
right-hand side of equ ations (I34bl) and (I34cl) are analogous to the corresponding terms in Cartesian 
coordinates derived in IgSOsI . and the remaining terms are geometric. 



4.5. Source Terms 



In summary, for the primitive variable equati ons in cylindrical coordinates, the MHD source 
term vectors are given by equations (18) and (19) of lGSOSi via the substitutions dxBx ^ R~^dR{RB}i) 
and dyBy ^ R~^d^B^^ and the geometric source term vector is given by 

1 /^,2 1 td2\ 

-^{v^VR- -pB^BR) 

. (35) 

-^'^cj.BR 
-^vrBz 

Since the geometric source terms arise directly from the scale factors in the i?-partials, we associate 
the geometric source term Sgeom exclusively with the i?-direction. Note that ||sgeom|l ^ in the 
limit of vanishing curvature, i.e. as i? ^ oc. 



^geom 



We emphasize that the geometric source terms in equation (f35l) are used only in obtaining the 
L/R states, not for the final FV update. Finally, the gravity source terms for the L/R states are 
given by the cylindrical coordinate components of — V$ in the momentum equation, and there is 
no gravity source term in the energy equation. 



5. Spatial Reconstruction 



In Athena, spatial reconstruction is p erformed in a directionally -split fashi on us ing; piecewise 
polyn omial app roxima tions as outlined in IColella fc Woodward! Il984l (hereafter ICWl ) , and IColella 
199d (hereafter IColellal). Here, we focus on piecewise linear and quadratic reconstructions, which 
yield second- and third-order approximations to smooth profiles, respectively. For a given coordi- 
nate direction, ^, we form the piecewise linear or quadratic reconstruction of each primitive variable, 
a(^), from the set {a^} of cell-centered volume-averages (including ghost-zones) at time t^, holding 
indices j and k fixed. In each case, we require for consistency that the volume-average of the 
reconstruction equal the volume-averaged data in the zth cell, i.e. 



(«({)). 



Vijk Jv.jk 



a(0 dV. 



(36) 



- 13 - 



Instead of defining a(^) in the ith zone explicitly, i.e. for ^ G [Ci-i/25 Ci+i/2]5 we find it more 
convenient to define the auxiliary parameter s E [0, 1] by 

s ^ (37) 

where = ^^+1/2 — is the width of the interval, so that ^ = Ci-1/2 + ^ 

We also employ slope-limiting and monotonization procedures to ensure that the resulting 
reconstructions are total-variation-diminishing (TVD) while providing somewhat steeper slopes at 
discontinuities. Of course, this can destroy the local formal order of the reconstruction, especially at 
extrema, but we pay this price for stability. Note, however, that while monotonicity is a sufficient 



condition for a r e const ruction to be TVD, it is not always necessary flLevequd I2OO2I ) . Recently, 



Colella fc Sekoral (120081 ) have described a slope-limiting method that, when combined with piecewise 
quadratic reconstruction, preserves the local order of convergence of the reconstruction at extrema. 
This has been implemented for Cartesian coordinates in the latest versions of Athena, but not for 
cylindrical coordinates, hence will not be described further here. 

For reconstructions in Cartesian coordinates, the procedures for the y- and z-directions are 
identical to the procedure for the x-direction. For the reconstructions in cylindrical coordinates, 
the only non-trivial difference from the Cartesian procedure occurs in the i?-direction since the 
discrete cell-volumes change with i?, but not with or z. Thus, we take ^ = x for the Cartesian 
cases and ^ = i? for the cylindrical cases. The Cartesian formulae apply, with suitable relabeling 
of coordinates, for 0- and ^-reconstructions. 

5.1. Piecewise Linear (2nd-order) Reconstruction 

Piecewise linear reconstruction approximates each primitive variable by defining in the ith 



zone, 



a{s) = aL,z + s Aai = an^i - (1 - 5) Aa^, (38) 



where Aa^ = aR^i — a^^i represents the difference of some quantity a over the zone, and a^^i and aR^i 
are the values of a at the left and right interfaces of the zone, respectively. Thus, to specify a{s) 
completely, we need only to define Aai and aL,i for each zone as functions of the volume-averages, 

5.1.1. PLM in Cartesian Coordinates 
From the consistency requirement in equation (136]) with Cartesian coordinates, 

(39) 



I rxi+i/2 

— / a{x) dx — I a{s) ds. 
^xJx,./o Jo 



Substituting equation (138]) into equation (139]) and integrating, we obtain 

di = ciL,i + ^Aai = aR^i - ^Atti, (40) 
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from which 



aL,i = cii- ^Aa^, (41a) 
dR.i = cii + l^cii- (41b) 



The usual Cartesian formulae for the differences over the zone are 



Aai = < 



I (a^+i — ai-i) , centered 

(a^+i — a^) , forward . (42) 
{ai — ai-i) , backward 



It is clear from equation (l42l) that constant and linear profiles are reconstructed exactly. In order 
to make the reconstruction TVD, these differences (or "slopes," as they are commonly called) are 
limited using; a slight variant of the monotonized central-difference (MC) limit er as described by 



Levequd (120021 ) . and flattened to avoid the introduction of new extrema. Leveque argues that lim- 
iting should be performed in characteristic variables so that the accuracy of the reconstruction for 
smooth wave families is not adversely affected by limiting in other non-smooth wave families. This 
requires a special description for systems of conservation laws including a bounded linear transfor- 
mation from primitive to characteristic variables, with inverse transformation from characteristic 
variables back to primitive. We do not go into detail here, but note that the inverse transformation 
is not guaranteed to be monotonicity-preserving, hence an additional monotonization is performed 
on the resulting primitive variable differences. In the previous implementation, this was done in 
a non-conservative manner, and we have since implemented a related scheme which preserves the 
consistency requirement of equation (|36l) . 

Finally, by applying the monotonized, limited differences Aa^ in equations (l40l) and (l4T1) . for 
smooth profiles we obtain second-order accurate approximations to the values of a across each zone 
at time t^, except possibly at local extrema. 

5.1.2. PLM in Cylindrical Coordinates 
From the consistency requirement of equation (f36l) with cylindrical coordinates, 

= r^^'^' <R)RdR C a{s)(R,_y2 + sAR) ds. (43) 

Substituting equation ([381) into equation (j43l) and integrating, we obtain 

ai = + \ Aai (1 + 7i) = aR^i - \ Aai (1 - 7^) , (44) 

where 

is a correction factor for curvature. Note that 7^ ^ for fixed Ai? as i? ^ 00, or for fixed 
R as Ai? 0, and from the formulae in equation (j44l) , we recover the Cartesian formulae in 
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equation (140]) . Solving for aL,i and aR^i, we obtain 

«L,i = cii- l^cii (1 + 7i) 1 (46a) 
= + I Atti (1 - 7i) . (46b) 

Next, we wish to define the differences Aa^ such that constant and hnear profiles are reconstructed 
exactly. Assuming a linear profile, say a{R) = CR with some constant slope C, we enforce the 
consistency relation of equation (143]) to obtain 

ai^{a),^{CR),^C{R),, (47) 

where 

{Rh = R^ + ^ (48) 

is the i?-coordinate of the volume centroid of the zth zone. To obtain an exact reconstruction 
for a linear profile, we require that Aa^ = C AR. First, we consider the Cartesian formula for a 
centered-difference to obtain 

i(a,+i - a,_i) = i (C(i?,+i) - C(i?,_i)) = CAR (l - ^^^^L^ . (49) 

If we divide the Cartesian centered-difference formula in equation (|42]) by the term in parentheses 
from equation (149]) . we obtain the desired difference, CAR, exactly. This — along with similar 
calculations for the forward- and backward-difference slopes — suggests the following definitions: 

I (tti+i - a^-i)/ (l - isi^fS-r) ' centered 

(a^+i - ai)/(l - i2r!1\ r, ) ' forward . (50) 

{ai — ai-i)/ f 1 — i2R^R. 1, backward 



Aai = < 



It is clear from equation ([501) that constant profiles yield Aai = 0, hence these are also reconstructed 
exactly. In order to make the reconstruction TVD, these differences are limited and monotonized 
as in the Cartesian case. Finally, by applying the resulting differences to equations (f38]) and (146]) . 
for smooth profiles we obtain second-order accurate approximations to the values of a across each 
zone at time t^, except possibly at local extrema. 

5.2. Piecewise Parabolic (3rd-order) Reconstruction 

Piecewise parabolic reconstruction approximates the profile of each primitive variable in the 
ith zone as 

a{s) = aL,i + s Aai + ^(1 - s)aQ^i = aR^i - {1 - s) Aai + s{l - s)aQ^i, (51) 



where, as for linear reconstruction, Aai = ciR,i — cll.i represents the average difference of some 
quantity a over the zone, and a^^i and aR^i are the values of a at the left and right interfaces of the 



-16- 



zone, respectively. The term a^^i is the so-cahed parabohc coefficient (see ICWi . eq. 1.5). To specify 
a{s) completely, we need only to define Aa^, aL,i, and ae,^ for each zone as functions of the set of 
volume-averages, {a^}. 

5.2.1. PPM Reconstruction in Cylindrical Geometry 

To satisfy the consistency requirement of equation (f551) in cylindrical coordinates, we substitute 
equation ([5T]) into equation and integrate, then solve for aQ^i in terms of the quantities and 

tt6,i = Q[ai- aL,i - ^Aa^ (1 + 7^)] . (52) 
E quation (|52|) is equ i valent to the defi nitio n of the parabolic coefficient appearing in equation (18) 



of iBlondin fc LufkinI (|l993l ) (hereafter IBLI ) . Recall that 7^ ^ for fixed AR as i? ^ oc or for fixed 
R as AR (see eq. HSj), and in these limits we recover the Cartesian version of the parabolic 
coefficient, which is the same as equation (l52]l except with 7^ = (see >CWi . eq. 1.5). 

We wish to define the values of aL,i, ciR,i, and Aa^ = aR^i — aL^i such that constant, linear, and 
parabolic profiles are reconstructed exactly, or a t least to second-order. By constructing a quartic 



polynomial from the {a^}, one can show fseelB N. ^3) that the cylindrical formulae can be obtained 



from the corresponding Cartesian formulae (see ICWi . eq. 1.6) by making the canonical substitution 
ai ai Ri. This yields 

aL,i Ri-1/2 = I (<^i Ri + Cii-l Ri-l) — I i^Cii Ri " Ri-l) -> (53a) 

aR^i i?i+i/2 = \ (<^i+l Ri^l + CLi Ri) — \ (5<^i+l Ri^l — ^CLi Ri) • (53b) 

Here, the centered-, forward-, and backward-differences in zone i are 

\ (a^+i - ai-i Ri-i) /Ri, centered 
6ai = <^ (a^+i Ri^i - aiRi) /i?i+i/2, forward , (54a) 
{ai Ri - ai-i Ri-i) /i?i-i/2, backward 

and similarly for the i + 1 and i — 1 zones. Note that the corresponding Cartesian formulae are 
given by taking i?^ ^ 1 in equations (|53]) and (|54|) . 

First, assuming a constant profile, a{R) — C, and taking the volume average across the zth 
zone, we have that ai — C. Using centered-differences, it follows that 6ai = C AR/Ri, hence from 
equations ([53]), we see that aL^i — C as desired. For forward- and backward-differences, a constant 
profile yields a^ = C [l + O {{ARf)] . 

Next, assuming a linear profile, a{R) — CR, and volume-averaging, we have given using 
equations (l47|l and (l48]l above. Using either centered-, forward-, or backward-differences, it follows 
that 5ai — 2C AR, hence from equations (155]) , we see that a^^i = CRi_i/2 and aji^i = CRi^i/2, as 
desired. 



Finally, assuming a parabohc profile, a{R) = CR^, it is straightforward to show that the 
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volume average of R across the ith zone is 

ai^C{R^).^C { Rf + 



my 



(55) 



Using centered-differences, it follows that Sai = C Ai? [Si?^ + 5(Ai?)^/(4i?^)] , hence from equa- 
tions (EHj), we see that aL,i = ^-^^-1/2^ desired. For forward- or backward-differences, it can be 

1 + (Ai?/i?^_i/2)^ • Thus, we conclude that parabolic profiles can 



shown that aL ^ 



^^1-1/2 



be recovered up to the required order. 

As in §5.H the slopes (meaning the centered-, forward- and backward-diff erence ScL j ^s) ar e 
monotonized in characteristic form using a slight variant of the MC limiter, as in Leveque fcood ) . 
Then, after computing the parabolic interpolant, th e slo pes are re-monotonized to ensure that the 
interpolation introduces no new extrema. Following [BIJ, the values of a^^i and aR^i are reset to 



whenever is a local extremum with respect to aL,i and aR^i, or to 

Qcii - ciR,i (4 + 37i) 



* 



2-37^ 
6ai - aL,i (4 - 87^ 



(56) 

(57a) 
(57b) 



2 + 37, 

whenever they are close enough to so that the parabolic interpolation function introduces new 
extrema. The test for this case, 

\ciR,i - CiL,i\ > |tt6,i|, (58) 

is geometry-independent. As a result, for smooth profiles we obtain second-order accurate approx- 
imations to the values of a across each zone at time t^, except p ossib ly at local extrema. By taking 
7i = 0, we recover the Cartesian versions of equations (f57l) (see ICWi . eq. 1.10). 



6. Characteristic Evolution and Averaging 

The final step in the calculation of the one-dime nsion al L /R state s is a characteristic time- 
evolution from to following the methods of ICWl and IColellal . This is accomplished by 
computing the time-averages of the solutions to the linearized primitive variable systems described 
in ^at zone interfaces over this half-timestep. The particular form of the averages depends on the 
direction, the order of the reconstruction, the coordinate system, etc. 



First, recall the modified primitive variable system described in §4] by equation (ETj), where A 
is the linearized hyperbolic wave matrix for the ID equations projected in the i?-direction, which 
is given by equation (123]) . Recall further that (weakly) hyperbolic systems of conservation laws 
have square wave matrices with real eigenvalues. Thus, let < • • • < represent the M real 
(but not necessarily distinct) eigenvalues of A corresponding to the M linearly independent left- 
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and right-eigenvectors, {/^,r^}, where v — 1,...,M. These eigenvectors are orthonormahzed so 
that • — Thus, any vector w G (in particular the centered-, forward-, or backward- 
differences across the zth zone, ^Wi) has the right-eigenvector expansion 

M 

^ = ^aV^ (59) 

with the coefficients oJ^ — - w representing the components of the projection of w onto the 
left-eigenspace of A. 

Next, we note that the characteristic form of the primitive variable system is obtained by 
multiplication of equation (j21l) on the left by L, the matrix whose rows are the left-eigenvectors 
of A, i.e. L = {/^, . . . ,/^}-^, hence LA = AL where A is the diagonal matrix consisting of the 
eigenvalues of A. Neglecting source terms for the moment, we obtain the homogeneous linear 
system 

dta + Kdna = 0, (60) 

where a = \jW is the vector of characteristic variables. From the form of equation (160 p . these eigen- 
values {A^} evidently represent the signal speeds of wave families along characteristics. Further- 
more, the system in equation (l60|) decouples into M constant-coefficient linear advection equations 
of the form 

dta + XOrq = 0, (61) 

which have the solution 

a(i?, r + r) = a^(C = i? - At), (62) 

where a^(^) is the reconstructed solution at time f^. Since the solution in equation (f62l) depends 
only on a(^), for each characteristic wave impinging on the interface, the contribution to the time- 
averaged interface state is given by the volume average of this reconstruction over the domain of 
dependence defined by the wave's characteristic speed. A, and the time interval (t^,t^ + At). 

In a time At, a right-moving wave travels a distance A At in the i?-direction. The volume 
in cylindrical coordinates of the domain of dependence of the left interface state at i?^+i/2 upon 
this wave is given by Fdod = (^i+i/2 ~ X At/2) X At Acj) Az. Thus, with XR = XAt/AR and 
XL = —XAt/AR^ we have the average of a over VboD equal to 

1 

fh+i/2(XR) = 7^ 1 ir^— / a{s) (i?,_i/2 + s AR) ds. (63a) 

Similarly, over the domain of dependence of the right interface state at Ri-1/2 upon a left-moving 
wave, we have 



1 f^^ 

fh-i/2(XL) = 7^ -1 / a{s) (i?,_i/2 + s AR) ds. 

' {Ri-l/2 + ^XL AR) XL Jo 



(63b) 
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6.1. PLM Evolution in Cylindrical Geometry 

Here, we describe the evaluation of the L/R states at time t^+^Z^ based on a piecewise hnear 
reconstruction of the underlying profile at time t^, defined by equations (l38ll . (l46ll and (l50ll of §5.11 
Substituting a{s) into equations (f551) and integrating, we obtain 

fL,i+i/2(XR) = - ^Xi? Atti (1 - /3i?,i(xi?)) (64a) 

on the right side of the zone (left of the interface) and 

fh-iMXL) = aL,^ + IxL Aa, (1 + (3lAxl)) (64b) 

at the left of the zone (right of the interface). Here, aL^i and aR^i are the values of a at the left 
and right interfaces of the zth zone, respectively, Aa^ is the monotonized difference of a across the 
zone from equation (|50]) , and we have defined the functions 

Mxr) = 77^ ^^^1^ An\ ^ 

6(i?i+l/2 - ^XR Ai?) 

PlAxl) - .^y (65b) 

6(i?i-l/2 + 2XL Ai?) 

as additional correction factors due to the curvature of the zone. Note that ^R^i{xR)^ h,i{XL) 
as i? ^ oo, i.e. in the limit of vanishing curvature, in which case equations ([Ml) reduce to the 
Cartesian formulae. Note further that in averaging a{s) over the whole zone, i.e. taking Xl/r — ^i 
we have PR^i{l) = /5L,i(l) = 7i (see eq. [45]), and from equations ([641) . we recover the averages of 
equation ([441) . 

6.2. PPM Evolution in Cylindrical Geometry 

Here, we describe the evaluation of the L/R states at time t^+^Z^ based on a piecewise parabolic 
reconstruction of the underlying profile at time t^, defined by equations ([5T1) . ([52l) . ([53l) . and ([541) 
of §5.2[ Substituting a{s) into equations ([63l) and integrating, we obtain expressions analogous to 
equations ([641) for the time- average of the right (left )-moving waves over the domain of dependence 
of the left (right) interface state at ^2+1/2(^2-1/2) upon these waves: 

/L,i+i/2(Xi?) = ^R,i - IXR [Aa^ - (1 - Ixr) Ci6,i] 

+ IXR [Atti - (1 - xr) I^rAxr)^ (66a) 
fh-i/2(XL) = aL,i + IxL [^ai + {1 - Ixl) ae^i] 

+ \XL [Aa^ + (1 - Xl) CLQ^i] PlAxl)^ • (66b) 



The functions ^R^i{xR) ^tnd ^L,i{XL) defined in equations ([65l) are the correction factors due to 
the curvature of the zone. The PLM result in equations ([64]) corresponds to setting ae,^ = 0. Note 
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that equations (l66l) also reduce to the Cartesian formulae (see ICWi . eq. 1.12) when (3 ^ as 
R ^ oo. Note further that in averaging a{s) over the whole zone, i.e. taking Xl/r — we have 
= /5L,i(l) = 7i5 and the results in equations are consistent with equation 

6.3. Sum Over Characteristics 



Once time-averaged L/R states have been obtained in the characteristic variables (as in eqs. 
[Ml or , we convert back to the primitive variables using R, the matrix consisting of the right- 
eigenvectors of A. Since L transforms from primitive to characteristic variables via a = hw, 
and since RL = I, where I is the standard identity matrix, Ra = w accomplishes the inverse 
transformation. 



Defining (^^£^^1/2 — fl i+i/2^^^B) ^ fRi-i/2(^L) each characteristic variable, 

we obtain the total time-averaged L/R states in primitive variable form by summing the projections 
of these contributions onto the right-eigenspace: 



w 



n+1/2 

L/R,i±l/2 



M 



n+1/2 

^L/R,i±l/2 ~ /_^'^L/R,i±l/2' 



(67) 



For example, in the Cartesian case fIStone et al.ll2008l . eq. 42), using equation fj64ap or fl66ap with 
;5 = 0, we have 



L,i+l/2 



At 
2A^ 



As writte n in equatio n the inverse transformation includes a sum over all waves. However, 
following ICWl and IColellal . contributions to it; on a given interface from waves that propagate in the 
opposite direction may be discarded to yield a more robust solution for strongly nonlinear problems, 
i.e. the waves are upwinded in the appropriate direction. Thus, for left-moving waves with A < 0, 
we may set Xl/r — ^ equation (|64al) to obtain ^^1/2(0) ~ ^R,i- Similarly for right-moving 
waves with A > 0, we may set Xl/r = in equation (|64bl) to obtain /^^_i/2(0) ~ ^L,i- 



Stone et al.l (j2008l ) have noted that this upwinding destroys the formal second-order conver- 



gence for smooth flows. In this case, the ID L/R states are accurate to the desired order if and only 
if we account for all waves, including those propagating toward the interfaces from the outside of 
the zone. With this approach, the sum in equation (l68l) includes all A^. This means that for a given 
zone we integrate an extrapolation of the local reconstruction profile over a domain of dependence 
that lies outside the zone. However, for non-smooth flows, we have found that this can lead to 
significant errors near discontinuities, since we may end up extrapolating the local reconstruction 
beyond a point of discontinuity into a region where it is no long er a p; ood approximation of the 
profile. Thus, for flows that may contain discontinuities, we follow ICWl and restrict to integration 
only over characteristics that propagate toward the interfaces from the interior of the zone, as 
previously described. However, we do not use the reference states for waves propagating away from 
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the interface. Thus, we compute the states using 

<7+i/2 = < + li^^\-^ E >^''r-i^^Ur'', (69a) 

K7-l/2-^^-l(^^)^-^ E Anr.(AH.K, (69b) 

for Cartesian PLM, which we refer to as "upwind-only" integration. For PLM in cyhndrical 
coordinates, factors (1 — /3i^,i(C^)) (1 + ^"^^ included in the sums in equations (I69ap 

and (169bp . respectively, using equations fl64ap and (164bp . Equations fl66ap and (166bp are used to 
obtain analogous expressions for PPM in cylindrical coordinates. 



7. Finite Volume Method 

The FV method uses approximate time- and area-averaged interface fluxes to update volume- 
averaged quantities. In cylindrical coordinates {R,(f),z), the differential volume element is 

dV = RdRd(t)dz, (70) 

and the flnite grid cell volume and interface area are: 

V^jk = R^ARA(|)Az, (71a) 

^R;i±l/2J,k = Ri±l/2^(P^Z, (71b) 

-4^; i,j±i/2,fc = ARAz, (71c) 

^z;^,J,k±l/2 = R^ARA(j). (71d) 

To derive the FV method, we integrate the system in equations ([2]) over the volume of a given grid 
cell, apply Gauss's Divergence Theorem, and integrate in time from to t^^^ to obtain 

^n+l _ ri^ ( p p 77.^+1/2 

/^n+1/2 _ ^n+1/2 



Ri Acj) V ^5 ^,j+l/2,fc 0; i,j-l/2,A: 

_ ^ 1^^+1/2 _ T^n+1/2 

i,j,A:+l/2 ^ i,j>-l/2 

+ At5;+^/^ (72) 

where Q represents the volume-averaged conserved quantities, F represents the time- and area- 
averaged fluxes, and S represents the time- and volume-averaged source terms. 

For cylindrical coordinates, in ^we rewrote the ^-momentum equation in a modifled angular- 
momentum preserving form in order to reduce the number of source terms on the right-hand side 
of the system. Therefore, when we apply the procedure in equation ([721) to equation (fT3l) , it yields 
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the FV update for the angular momentum, Rpv^, not the hnear momentum, pv^. However, it can 
be shown that the radial contribution from a quasi-FV update of the pv^ equation (fT4l) . 



( p2 pn+1/2 _ p2 pn+1/2 \ .„„s 



is equivalent to the corresponding terms in the true FV update of equation (lllbp . 



At 



i?, Ai? 



/ pn+l/2 p pn+1/2 \ a+ZME^Y 

[Ri+l/2-t'R. i+l/2J,k - ^i-l/2J^R; i-l/2,j,k ) \~R~ / ■ k 



n+1/2 



(74) 



to second-order away from the origin for smooth flows. Note that equation (l74|) contains the volume- 
and time-averaged geometric source term, —Mji^f^/R — —(pvRV(f) — BjiB^f))/ whereas equation (173]) 
has no source term. 

The only nonzero component of the geometric source term, S^^^'^-jj^, required for the FV 
update is in the radial momentum equation flllap . For this term, we must compute 

Naively, one might compute this source term from volume averaged quantities at time f^. However, 
to achieve second-order accuracy, it is necessary to use time-centered estimates of these quantities, 
i.e. advanced to the half-timestep One can think of this as a sort of trapezoid rule applied 

to the time domain. This half-timestep advance is performed using a combination of FV updates 
on the volume-centered variables p and pv^j) and CT updates on the interface-centered B^j^. For the 
P* contribution, the FV and CT updates to t^+^Z^ are too costly, since they would be required for 
every variable and they, in turn, would require source term calculations. Instead, we compute the 
total pressure contribution directly from the fluxes at i?-int erf aces. The appropriate second-order 
average is 

' P* \ _ ^i+l/2^*+l/2 + ^i- l/2-F^-i/2 ^^g^ 



^ / ijk 

where time-averaged pressure fluxes returned directly from the Riemann solver. 

Our application of the geometric source term is similar to what is done for the gravity source 
terms in the existing Athena code. However, we have found it easier to maintain centrifugal 
balance numerically by using the analytic gravitational acceleration function g{x) = — V$(cc) in 
the momentum equation, rather than approximations of the gradient using finite-differences of the 
static potential, We approximate the gravitational source term for the momentum equation by 

SZl%k - -(PV$)S'/^ {P)l'k"^'9i{xh,,), (77) 

where {x)ijk is the volume-centroid of cell {i,j,k), the radial component of which is given by 
equation (l48l) . Note that for the case of solid-body rotation with uniform density, g (x so that 
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the gravitational source term given by equation (1771) is exact. For the energy equation, we rely 
on the previously implemented FV update based on the potential function which allows the 

energy equation to be written conservatively. 

To compute the gravitational source terms appearing in the calculation of the L/R states, 
which only appear in the momentum equation in primitive variable form, we evaluate g{x) at 
the area-centroid of each interface. Note that the i?-coordinate of the area-centroid of (/)- and 
^-interfaces coincides with the i?-coordinate of the volume-centroid of each adjacent grid cell. 

8. Constrained Transport 

In this section, we discuss m odifica tion s in cylindrical coordina tes to the constrained trans- 
port (CT) algorithm described in IgSOsI and Evans fc Hawley ( 1988 ). As argued in those papers, 
the integral form of the induction equation (lldjl is most naturally expressed in terms of finite 
area-averages rather than volume-averages. In this way, the equation becomes a statement of the 
conservation of total magnetic flux through a given grid cell and as such automatically preserves 
the V • B constraint. 

8.1. Integral Form and Consistency Relations 

To see this, we rewrite the induction equation as 

dtB + V X £ = 0, (78) 

where £ = —v x B is the electric field in ideal MHD (the electromotive force [EMF]). Then, 
integrating over the oriented bounding surface of grid cell {i,j,k) and applying Stokes' Theorem, 
we find that 

TDTi+l _ on fcn+l/'2 _^n+l/2 

^R- i±l/2,j,/c - ^R] i±l/2,j,k R.j^^i^ V ^' ^±l/2,j+l/2,fc ^z- i±l/2,j-l/2,/c 



At/ ^+1/2 _^n+l/2 \ . X 

Az V ^±l/2,j,A:+l/2 ^0; i±l/2,j,it-l/2y \i^c%) 



^0; ij±l/2,k ~ ^(/>; i,j±l/2,k ^ \^z- i+l/2,j±l/2,/c i-l/2j±l/2.k ) 

_ At / n+1/2 _ ^n+1/2 \ .^Q. X 

\^R] i,j±l/2,/c+l/2 ^R- i,j±l/2,/c-l/2 y ' yi^v^) 



TDU+l - m - ( pn+1/2 _ ^n+1/2 \ 

^z]i,j,k±l/2 ^z;iJ,k±l/2 R./\R \ (t>-.i+l/2,j,k±l/2 ^ cf)- i-l/2j,k±l/2 ) 

Ri A(f) V hj+i/2,k±i/2 ^R- i,j-i/2,fc±i/2 y ' y^^^) 
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where 

Bl;i±i/2,j,k ^ ^ X-TT- / / BR{Ri±y„cP,z,e)Ri±i/2d<Pdz, (80a) 

B;,^,J±l/2,k ^ XWX; / B^(R,ct>^^y,,Z,ndRdz, (80b) 

LAnLAZ Jzk_^/2 JRi-1/2 

I rRi+1/2 

Bl^,J,k±l/2 ^ B ABA6 / / B,{R,<j,,Zk±i/2,nRdRd<j,, (80c) 

i J (t)j-i/2 ^ Ri-i/2 

are the interface area-averaged components of the magnetic field normal to each surface (i.e. the 
magnetic flux per unit area) and 



^R;ij±i/2,k±i/2 = AtXR £R{R,(t)j±i/2,Zk±i/2,t)dRdt, (81a) 

^J^±1Aj>±1/2 = Ati?,^i/2 A(/) / £:^(i?,±i/2,0,2:fc±l/2,^)^z±l/2rf0rf^, (81b) 

+1/2 1 f^^^^ 

^."^iiAjiiAit = y^^ ^ f,(i?,±i/2,(/>,±i/2,^,^)d^dt (81c) 

are the corner-centered EMFs averaged over the edges bounding each surface. The EMFs in 
equations (l8T1) are approximated to some desired order of accuracy and the surface-averaged field 
components are evolved using equations ([79]) . 

The interface-centered, area-averaged magnetic field components in equations (fSOl) comprise 
the fundamental representation of the magnetic field in Athena. However, one often needs to refer 
to the cell-centered, volume-averaged magnetic field components as well. Therefore, we adopt the 
averages 



^z,j,/c = \^i-l/2^R- i-l/2,j,k ^ ^i+l/2^R- i+1/2. 



B^fh = TTTT ( /2^R. .--I /9 ^ + /2^r. .-^1 /2,j,/c) ' (82a) 



^tj!k = 2 (^^5 ^,j-l/2,fc + ij + l/2,k) ' (82b) 

Note the use of an i?- weighted average of the Br interface values in equation (|82al) ; it is straightfor- 
ward to show that this is the appropriate second-order accurate average in cylindrical coordinates. 
Equations (182]) imply consistency relations between the Godunov fluxes computed by the Riemann 
solver (the fluxes of the volume-averaged magnetic field components) and the corner-centered EMFs 
(the fluxes of the area-averaged rnagnetic field components). These relations define how they are 
computed from each other (see IGSOSL for details). 

8.2. Calculating the EMFs 



The primary modification to the CTU+CT algorithm described in iGSOSi for cylindrical coor- 
dinates concerns the calculation of the upwinded, corner-centered EMF component Sz- As we shall 
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demonstrate, we must combine spatial gradients of different curvature to form £z- However, to form 
£r or £(f)^ we combine spatial gradients of the same curvature, hence the effect of that curvature is 
cancelled out and no subsequent modification is necessary. 

For example, to compute £z;i-i/2,j-i/2^ we estimate (5(/)£^z)i-i/2,j-3/4 ^^^d use a centered- 
difference scheme to calculate one estimate: 



i-l/2,j-l/2 - £z- i-l/2,j-l n 2 V^0^2;;^_i/2j_3/4 ' 

In the same manner, we integrate £z to the corner from each of the remaining adjacent interface 
centers and take the arithmetic average: 

^z] i-l/2,j-l/2 = \ {£z; i-l/2J-l + ^z; i-l/2J + ^z; i-l,j-l/2 + ^z; ij-l/2) 

R^-l/2 M 



+ 



{d^£z 



(83) 



+ - 



+ 



{dR£ 



(^(/)^^2:)^_l/2j-3/4 (^</>^^)i-l/2,j-l/4 
2;)7:-3/4,j-l/2 ~ (^^^^)i-l/4,j-l/2 



(84) 



To ensure stability, for (5(/)<?z)i-i/2,j-3/45 we use the upwinding scheme (|GS05l , eq. 50) based 
on the sign of the mass flux at the center of each interface: 



{d(t)£z)i-l^j-?,/Ai 
{d(t)£z)i-l/2,j-?,/A = { {dcf)£z)ij-3/4, 

\ [(^(/>^^)i-l,j-3/4 + (^(/)^^^)i,j-3/4] 

The formulae for the remaining gradients are analogous. 



^i?; i-l/2,j-l > 
^R- i-l/2,j-l < 

otherwise 



(85) 



To obtain the estimates of d(i)£z needed in equation ([85l) , we use a centered-difference scheme 
based on the cell-centered EMFs, computed using volume-averages of p, pv^ and B, and on the 
interface-centered EMFs, which come directly from the fluxes: 



(^</>^^)i,j-3/4 ~ A(/) (^2:; i,j-l/2 £z] 1,3-1) 5 



(86a) 
(86b) 



Note that this scheme has no dependence on the particular type of Riemann solver used to calculate 
the fluxes. Furthermore, note the different radial scale factors appearing in equations (l86l) resulting 
from the combination of ^-gradients at different radii; the factors of A(/) cancel the factor of A(/) 
from equation (l84l) . but the radial scale factors themselves do not cancel and must be inserted into 
the algorithm. On the other hand, when ^-gradients are combined at the same radius, as is the 
case for the corner-integration of £r^ both the and the corresponding radial scale factors cancel, 
hence no modification is required. 



9. The Athena Algorithm 



In this section, we summarize in some what greater detail the main steps of the six-solve version 
of the CTU+CT algorithm adapted from IStone et al.l (120081 ) for cylindrical coordinates (see also 
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GS05I : IGSOSI . for details). 



1. Compute the first-order source terms, S^j j^, in both conservative and primitive variable forms 
using the initial volume-averaged data at time f^. These include the geometric source terms, 
'S'geom (eq. [35] for primitive variables, and eq. [75] for conserved variables), the gravitational 
source terms, S'grav? computed from static accelerations and potentials (eq. [77]for conserved 
variables), and the MHD source terms arising from the V • B constraint (e.g. eqs. [32] and [331 
for 3D). 

2. Compute the L/R interface states, Qf/_f/2jk^ using the desired 
reconstruction scheme on the initial data in primitive variable form. This requires reconstruc- 
tion with characteristic evolution as described in §^ and [6l followed by application of the 
parallel components of the (primitive variable) source terms from step ([T]). 

3. Compute the first-order interface fluxes, i^^. ^_i/2 j /c' i j-i/2 -^z- i j k-i/2^ from the 
interface states via an exact or approximate Riemann solver. 

4. Compute the corner-centered electric field components, i j_i/2 k-i/2^ i-i/2 j k-i/2^ 

^z' i-i/2 j-i/2k^ from components of the interface-centered fluxes from step (|3]) and the cell- 
centered electric field computed using the initial data at time t^, via equations (l84|) and (185]) 
for the z-components. 

5. Update the interface magnetic field components for a half-timestep using the CT difference 
equations ([79]) and the EMFs from step 

6. Compute the updated L/R interface states, Q^Lij2^]l^ ^ Q^^j^'1/2)!^ ^ 

applying transverse flux gradients to the non-magnetic variables of the interface states and 
then adding the transverse components of the source terms from step ([T]). 

7. Use the fluxes from step (|3]) and the source terms from step ([T]) to compute the velocities at 
the half-timestep using conservative FV updates of the cell-centered density and momentum 
at time t^. Average the half-timestep interface magnetic field components from step ([5]) to 
obtain the cell-centered magnetic field components at time t^+^Z^ using equations (l82]) . Then, 
calculate the cell-centered electric field components, £^l^j\^ ^z-\^jk^ using the 
cell-centered velocities and magnetic fields at time 

8. Compute the second-order interface fluxes, j k^ ^^^/j-i/2fc' -^^^lj'k-1/2^ using 
the updated interface states from steps ([5]) and ([6]) via the Riemann solver. 

9. Compute the corner-centered electric field components, /c-i/2' ^^•^i^i/2 j k-1/2^ 
^z-i-i/2 j-1/2 k^ from components of the updated fluxes from step (|H]) and the cell-centered 
electric field components computed in step ([7]), as in step (|4]). 
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10. Use the fluxes from step ([Hj) to obtain half-timestep conservative FV updates of the ceh- 
centered density and ^-momentum at time t^, similar to step ([7]). Compute the cell-centered 
total pressure, P*, from the interface-centered quantities returned by the Riemann solver in 
step (jSj) using equation (I7B]) . Combine the cell-centered B^f) from step ((Tj) with p, pv^f), and 
P* to construct the second-order geometric source term (eq. [75l) at time Combine 
p and the static gravitational acceleration, gr, to construct the components of the gravity 
source term for the momentum equation at time t^^^/^. Combine the interface-centered pv, 
obtained directly from the fluxes in step (|5j), to construct the gravity source term for the 
energy equation at time 

11. Using the fluxes from step (|8j) and source terms from step (ITO]) . advance the cell-centered 
quantities from time to t^^^ using conservative FV updates on the hydrodynamic variables 
(mass, momentum, and energy) and using the CT difference equations (l78]) with the EMFs 
from step (|9]) to update the interface magnetic field components. 

12. Average the updated interface magnetic field components from step (|TT]) to compute the 
updated cell-centered values using equations (182]) . 

13. Increment the time to t^^^ = + At and then compute a new timestep using the standard 
CFL condition based on the maximum signal speed at cell centers and on the size of the grid 
cells. Here, we must use Ri Ac/) to estimate the CFL stability criterion, since the Riemann 
solvers compute linear wavespeeds. 

This algorithm is simplified for the purely hydrodynamic case. Besides having fewer variables 
to store, reconstruct, and evolve, there is no need to compute MHD source terms, magnetic compo- 
nents of geometric source terms, or corner- or cell-centered EMFs, or to apply FV or CT updates 
to magnetic field components. 

In adapting the code for cylindrical coordinates, we have altered several steps of the original 
Cartesian algorithm to varying degrees. First, we significantly change the computation of the L/R 
states for the P-direction in step ([2]); only a minor change is needed for the ^-direction and no 
change is needed for the z-direction. Second, we add geometric scale factors to the flux differences 
in the conservative FV updates performed in steps ([7]), (|T0]) . and (|TT|) . and we include the geometric 
source terms computed in steps ([T]) and (|T0]) . As detailed in §^and[4l the geometric source terms 
applied in steps ([2]) and (|6j) differ from those applied in steps ([7]), (|T0]) . and (fTT|) because of the 
differences in the primitive and conservative forms of the evolution equations, and furthermore, the 
source term contribution from the total pressure, P*, is computed differently in steps (P) and (|T0]) . 
as discussed in ^ Finally, we change the CT calculation in steps (j4]) and (ITO]l to reflect the 
additional geometric scale factors appearing in the cylindrical coordinate version of the induction 
equation (|78]) and to enforce modified consistency relations, as described in §8.11 
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10. Code Verification Tests 



In this section, we present a suite of tests of our cylindrica l coordinate adaptation of the Athena 
code . Some are dra wn from tests pubhshed by other authors (I gSOS . GS08 . Londrillo fc Del Zanna 
2OOOI . ISakurailll985l ). which were originally written to test Cartesian codes, including Athena itself, 
while others are new. Where possible, we have tried to make comparisons with the existing Carte- 
sian tests in order to demonstrate our code's ability to recover their results both qualitatively and 
quantitatively. We include tests in one, two, and three spatial dimensions, in both hydrodynam- 
ics and MHD, with solutions that are both smooth and non-smooth, and having varying levels of 
symmetry. 



10.1. Force Balance 



In this problem, we investigate various steady equilibria in order to evaluate the code's ability to 
balance forces. While there are many possible tests to choose from, we give here two representative 
examples demonstrating simple magnetohydrostatic equilibria. 

First, we consider the axisymmetric magnetic field 

B = (87) 

for which the outward magnetic pressure and inward tension forces sum to zero. Note that the 
magnetic field given by equation (f87l) satisfies the V • S = constraint. We use Bq = 1 and set the 
velocity 1; = 0, mass density p = 1, and gas pressure P = 1. Figure [T] shows the convergence of the 
L2 norm of the Li error vector (RMS error) for the solution at time t = 10, defined as 

i 

where is the initial solution. For reference, we plot a line of slope —2 (dashed) alongside the error 
to demonstrate that the convergence is second-order in These data were computed using the 

HLLD fluxes and third-order reconstruction; the results were similar for all combinations of Roe 
or HLLD fluxes, second- or third-order reconstruction, and ID, 2D, or 3D integrators. However, 
because the 2D and 3D algorithms differ significantly from the ID version, especially in their 
inclusion of transverse flux gradients and CT updates, we also present the results of the same test 
using these integrators on grids which are essentially one-dimensional, but contain a few grid cells 
in each transverse direction considered. By symmetry, it is clear that any number of cells may be 
used in the transverse directions, but since the grid cell volumes change with R in multidimensions, 
the CFL condition will determine the timestep based on grid cell volume as well as the maximum 
signal speed, so the absolute errors should not be compared between, say, the ID and 2D algorithms. 
Only the order of convergence of each individual algorithm is meaningful. Additionally, we have 
performed tests with the outward acceleration v'^/R from a solid-body rotation profile V(f) = QqR 
balanced by the gradient of the static gravitational potential $ = (fio^)^/2, as well as with constant 
Vz 0, and find similar results in all cases. 
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Second, we consider the non-axisymmetric magnetic field 



B = ?^^R^ (89) 
which when combined with the gas pressure 

yields zero net force in the ^-direction. The combination of static gravitational potential 
and mass density 

p = Po[l + sin2(V;)] (92) 

together balance the gradient of the gas pressure in the i?-direction. Here, we use the angular 
coordinate '0 = 27r((/) — 0min)/(0max — 0min) so that Br^ P, and p are all periodic in the ^-domain. 
Note that the magnetic field given by equation f[5^ satisfies the V • S = constraint. We use 
So = 1, ^0 = 1, Po = 1, and once again use the solid-body rotation profile v^f) — QqR balanced by 
the gradient of the static gravitational potential $2 = (f^o^)^/2, where Qq = 7r/4. Note that the 
total potential is given by $ = $1 + $2- 

Figure [2] shows the convergence of the RMS error for the solution at time t = 10. For ref- 
erence, we once again plot a line of slope —2 (dashed) alongside the error to demonstrate that 
the convergence is second-order. These data were computed using the HLLD fluxes, third-order 
reconstruction, and the 2D integrator; the results were similar for all combinations of the Roe or 
HLLD fluxes, second- or third-order reconstruction, and 2D or 3D integrators. We use a computa- 
tional domain of size N-hj-N with R G [1,2], (j) G [0,7r/4], and 2: = 0. We use periodic boundary 
conditions in the (/)- and z-directions, and Dirichlet boundary conditions in the i?-direction. Similar 
results were obtained using a Neumann boundary condition in the i?-direction. 

10.2. Rotational Stability 

In this problem, we investigate the stability of rotating disks evolved with our code, using 
the 2D integrator. Given a differential rotation profile, Rayleigh's criterion for stability to 

axisymmetric, infinitesimal disturbances is that specific angular momentum increase outward: 

dR [{R^n{R)f] > 0. (93) 

While it is possible that systems satisfying Rayleigh's criterion are still subject to growth of finite- 
amplitude non-axisymmetric disturbances, laboratory measurements at Reynolds number up to 
2 X 10^ have found that Couette flows violating Rayleigh's criterion show large angular momentum 



trans port associated with turbulence, while those satisfying Rayleigh's criterion do not (jji et al 



20061). 
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Fig. 1. — Convergence of the RMS error in the Li-norm for the force-balance problem in ID, 
2D and 3D. For reference, we have plotted a line of slope —2 (dashed) to show that the convergence 
is second-order in 




Fig. 2. — Convergence of the Li-error for the 2D Br force-balance problem in 2D. For reference, 
we have plotted a line of slope —2 (dashed) to show that the convergence is second-order in 
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For our test, we consider power-law rotational profiles of the form Q{R) oc where g is a 
constant, the so-called "shear parameter." Rayleigh's criterion in a differentially rotating system of 
this form predicts stability when ^ < 2, and instability for q > 2. For example, Keplerian rotation 
with g = 1.5 is predicted to be stable (for unmagnetized flows). 

Using a constant background density and pressure, we set 

v^(R) = Rn{R) = noR^-"^ (94) 

and set the gravitational potential so that rotational equilibrium is achieved. Next, we perturb v^j) 
to 

Vcf) = Vcf> + Svcf>, (95) 

where Svcf) is a random variable uniformly distributed in [— e,e], and e is small, typically on the 
order of 10~^. The same initial perturbation is used for each value of q considered. We use a grid 
of 200 X 400 cells over the domain [3,7] x [0,7r/2], which is chosen so that i?avg ^ 2AR. We 
set po = 200, Pq = I5 and use an adiabatic index of 7 = 5/3 which gives Cg ^ 0.09. With = 27v, 
'^(/),min ^ 0.9, which puts the Mach numbers in a range of approximately 10-20 over the domain. 
Thus, the flow is rotationally dominated. 

As a diagnostic of instability, we compute a scaled mean perturbed angular momentum flux, 

{RpvR Sv^) _ J J RpvRivd, - RVt)RdRd(j) 

{RP) ~ JjRPRdRdct) ' ^ ^ 

For stable flows, this will remain on the order of the initial perturbation, but for unstable flows, it 
will diverge exponentially. Figure [3] shows the values of the dimensionless angular momentum flux 
as a function of time for t G [0, 300] for various values of the shear parameter near the marginal 
stability limit of g = 2. Consistent with Rayleigh's criterion, the flows with q < 2 remain stable, 
and those with q > 2 go unstable. For the q = 2.05 case, the instability reaches saturation more 
quickly (around t = 90) and the mass flies off the grid, but the characteristic exponential growth 
is observed before this point. This test demonstrates the code's accurate conservation of angular 
momentum near the boundary of rotational stability. 

Additionally, we investigate the long-term stability of the rotation profiles given by the shear 
parameters q — 1^ typical of galactic disk systems, and q — 1.5, typical of Keplerian systems. 
For this test we use the unperturbed equilibrium solutions as initial data. As a diagnostic of the 
error, we compute the cumulative mean of the dimensionless background angular momentum flux 
(proportional to the radial accretion rate): 

{{RpvRiytR))) _ JJJ RpvR{VtR)RdRd(i)dt 
{{RP)) ~ /// RPRdRd(j)dt 

Figure m shows the dimensionless cumulative mean as a function of time for t G [0, 300]. We use the 
same computational domain and background state as above. Since V(j) — Q.R oc R^~^ is constant 



(97) 
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for g = 1, and constant profiles are reconstructed exactly in our code, we observe relatively small 
errors for this level of discretization, indicating that angular momentum is conserved very well for 
systems of astrophysical interest. 

10.3. Adiabatic Blast Wave 

I n this problem, we investigate a strong 2D shock using the HLLC solver. We use the parameter 



set of IgsosI and compare the outputs of our cylindrical code with those from the Cartesian version 



of Athena. For the Cartesian version we use the domain G [—0.5,0.5] x [—0.75,0.75], and 

for the cylindrical version we use the domain (i?, 0) G [1,2] x [—0.5,0.5] so that the physical 
domain spans an arc-length of i?mid(0max — 0min) = 1-5 at i?mid — 1-5, giving a roughly similar 
domain sizes. The initial conditions consist of a circular region of hot gas with radius i? = 0.1 
and pressure P = 10 in an ambient medium of uniform pressure Pq = 0.1 and density po = 1- 
We use a computational grid of 200 x 300 cells, third-order reconstruction, and the HLLC fluxes 
with upwind-only integration for the L/R states in each version of the test. Contour plots of 
the density, pressure and specific kinetic energy densities of the evolved state at time t = 0.2 are 
shown in Figure [5] using the cylindrical (left column) or Cartesian (right column) version of Athena. 
For additional comparison, ID plots of these variables along a horizontal line through the center 
of the blast are shown in Figure [6l demonstrating excellent agreement between the Cartesian and 
cylindrical versions of Athena. Notice in the Cartesian version that symmetry is perfectly preserved 
by the integrator, which is most easily seen in the grid noise in the interior of the shell in Figure El 
Symmetry is also preserved rather well in the cylindrical version although in this case the grid is 
non-uniform. 

10.4. Rotating Wind 

In this problem, we investigate a steady, axisymmetric, rotating hydrodynamic wind as a test 
of angular momentum transport across a sonic transition. We adopt a Newtonian gravitational 
potential of the form = —GM/R. The constants of motion are given by 

K = Pp-\ (98) 
M = RpvR, (99) 
J = Rv4,. (100) 

This flow must satisfy the BernouUi equation, B — constant, along streamhnes for 

B = \v'^ + h + ^g, (101) 



where 



is the specific enthalpy of the gas. 



h^l^ = ^ (102) 
P 7-1 



We scale density and pressure to their values at infinity, /?oo = p(oo) and Pqo = -P(oo)) and 
the radial coordinate to some finite fiducial value, Rb- In terms of a = p/poo and x = R/Rb, the 
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Fig. 3. — Mean dimensionless angular momentum transport as a function of time in the Rayleigh 
rotational stability test for various values of q. 
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Fig. 4. — The cumulative mean (time- and space-average) of the dimensionless background radial 
angular momentum flux as a function of time for shear parameters q = 1 and g = 1.5 with 
unperturbed initial data. 
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Fig. 5. — Contours of selected variables of the evolved state (at time t = 0.2) for the 2D hydro- 
dynamic blast wave test using 200 x 300 grid cells, third-order reconstruction, HLLC fluxes, and 
the cylindrical (left column) or Cartesian (right column) versions of Athena. Thirty equally spaced 
contours between the minimum and maximum are drawn in each plot. 
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Fig. 6. — Plots of selected variables along a horizontal line through the center of the blast at time 
t = 0.2 for the 2D hydrodynamic blast wave test (see Fig. [5] legend), using the cylindrical (circles) 
or Cartesian (solid line) versions of Athena. 
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constant entropy parameter is K = c^/(7Poo where = j Poo /poo is the square of the sound 
speed at infinity. At any radius, the local sound speed, radial velocity, and specific enthalpy satisfy 

c2 = c^a^-\ (103) 
Vr = Mlcla^-\ (104) 
h = -l^c^a^-i, (105) 

where Mr = vr/cs is the radial Mach number. 

We define the dimensionless radial mass flux by 

A ^ ^ ^ = xMna'^'^+'y' (106) 

-t^BPooCoo 

and the dimensionless angular momentum by 

^B^oo ^oo 

Solving equation (11061) for a and introducing (3 = 2(7 — l)/(7+l), we find that a^~^ = [^/{X-^r)]^ - 
Finally, taking Rb to be the Bondi radius, Rb = GMjc?^^ we obtain the dimensionless Bernoulh 
equation 



—y + xO-^-iry!'-'' 

7 — 1 1 



(108) 



where ;B = B/hoo is the dimensionless Bernoulli constant, and hoo = c^/(7 ~ 1) is the specific 
enthalpy at infinity. 

Figure El shows the contours of A for various values of cj, using an adiabatic index of 7 = 5/3. 
For the ^ — ^ ca se, we recover a cylindrical version of Parker's spherically symmetric wind (see, 
e.g., Sit zerlll978l ). The bold lines represent the transonic solutions (wind and accretion) passing 



through the X-type saddle points. These solutions can be found by first writing equation (11081) as 
J-'{M.R)f(X) = G{x) requiring that J^' — Q' — 0. The first constraint implies that Mr — 1, 
i.e. the saddle point is the sonic point. The second constraint yields a quadratic in and for 
uj G (0,a;max)5 there are two distinct, positive solutions. 



3-7±A/(3-7)2-16i3cj2 
X± = ^-^--^ ' (109) 

where c^max = (3 — 7)/ (40-'^/^). It can be shown that X- gives an 0-type critical point and x+ gives 
the desired transonic critical point. For uj > c^max, no transonic solutions exist. Note further that 
no transonic solutions exist for x < Xmin, where Xmin represents the point at which the transonic 
wind solution for which Mr ^ 00 diS x ^ ^ joins the transonic accretion solution for which 



-37- 



A4r ^ as X ^ CO (see, e.g., the lower-left panel of Figure [7j). We define Xmin to be the smallest 
value of X ^ for which G{x) ^ 0? which is given by 



-(7-1) + V(7- 1)2 + 2^^^2 

Xmin = . (110) 



Finally, the critical value of the radial mass flux, Ac, is defined by equation (llOSp with x — X+ 
yWi^ = 1, which is given by 



xi-\x^-^'f'. (Ill) 



For our code test, we use 7 = 5/3 (i.e. /3 = 1/2), cj = 0.3, and 5=1, which give the transonic 
solution shown in Figure Eb. The critical point occurs at x+ ^ 0.479, with Ac ^ 1.377. We solve 
the problem on the domain x ^ [X-?^], where X- ^ 0.188, using bisection with a tolerance of 
e = 10~^^ to evaluate Mr at each x from equation fllOSp . Once Mr is known, vr, p, and P follow 
algebraically. We choose units such that GM = Coo = 1, which yields Rb = I- We fix the solution 
at the inner and outer boundaries and evolve the initial solution long enough for it to settle into 
equilibrium. 

Figure [8] shows the convergence of the L2 norm of the Li error vector for the solution at time 
t = 5.0. These data were computed using the Roe fluxes, second-order reconstruction, and the 
ID; the results were similar for all combinations of Roe or HLLC fluxes and second- or third-order 
reconstruction. The test was also performed using the 2D and 3D integrator with a few grid cells 
in the transverse directions. 

This test clearly demonstrates the code's ability to maintain smooth, steady hydrodynamic 
flows in both subsonic and supersonic regimes, as well as its ability to conserve angular momentum 
to second-order in cylindrical geometry. 

10.5. Field Loop Advection 

In this problem, we investigate the advection of a w eak fie l d loop in 2D and 3D cylindrical 



coordinates, analogous to the Cartesian test appearing in IGS05I : IGS08I . The main difference with 
our test is that we advect the field loop in the ^-direction only as opposed to a more general 
advection oblique to the grid. We use the computational domain (R, (j)) G [1, 2] x [—2/3, 2/3], which 
has the same total area as the Cartesian version of the test. We use periodic boundary conditions 
in (j) and fixed boundary conditions in R. We use uniform initial density po = 1 and pressure Pq = 1 
with a solid-body rotation profile of v^i^ = fio^, where we set = 4/3 so that the field loop is 
advected once across the grid hj t — 1. The initial z-component of the magnetic field is 0, and the 
R- and ^-components are set using the z-component of the magnetic vector potential 

f ^o(ao-r) forr<ao, 

I for r > ao ' ^ ' 



-38- 




Fig. 7. — Contours of the dimensionless mass flux, A, for a rotating hydrodynamic steady flow with 
7 = 5/3 and dimensionless Bernoulli constant B — 1. The scaled angular momentum in each panel 
is (a) = 0, (b) = 0.2, (c) uo = 0.3, and (d) = 1/3. The critical transonic contours are shown 
in bold. 
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where ag is the radius of the field loop and r = ^i?^ + i?Q — 2RRq cos((/) — ^o) is the distance from 
the center of the loop (i^o^^o)- We use Aq = 10~^ and ao = 0.3 so that inside the field loop 
P/B = 10^ and the field loop should be advected passively. 

Figure [9] shows the magnetic energy density and magnetic field lines at times t = and 

t — 2 for the 2D problem. The field lines are the contours of Az which can be readily computed 
since the field is planar and the CTU+CT algorithm preserves the V • S = condition. Note that 
the circular shape of the field lines is nicely preserved. 

Figure [10] shows the time evolution of the volume-averaged magnetic energy density. The 
dissipation is well-described by a power law of the form {B^ /2) = ^4(1 — (t/r)^) (for t <C r), with 
A = 7.02 X 10"^, T = 1.46 X 10^, and a = 0.342, and with a residual error of 0.0580. Note that the 
overall dissipation in this pr oblem is less than that of the Cartesian version since the advection is 



only in one direction; IGSOSI found r = 1.06 x 10^ and a = 0.291 for a similar fit. 



10.6. Blast Wave in a Strong Magnetic Field 

In this problem, we investigate 2D and 3D MHD shocks in a strongly magnetized medium 
with low plasma-/?, denote d 0p = 2P/B. We run two problems, one with Bq = 1 and Pp = 0.2 
using the parameter set of IGSOSI. a nd another with = 10 and /3p = 0.02 using the parameter 



set of iLondrillo fc Del Zannal (120001 ) . In each case, we compare the outputs of the cylindrical and 



Cartesian versions of Athena. 

The moderate S-field, /3p = 0.2 case uses HLLD fluxes and the same setup as the hydrodynamic 
blast described in §10. 3| but with a uniform background magnetic field of strength = 1 oriented 
at a 45° angle to the positive x- or i2-axis. For the 2D case, contour plots of the density, pressure, 
specific kinetic energy, and magnetic energy are shown in Figure [TT] based on the cylindrical (left 
column) or Cartesian (right column) versions of Athena. Figure [T2l shows a plot of these variables 
along a horizontal line through the center of the blast. Note that in the Cartesian version, the 
background field is uniformly inclined to the grid, but in the cylindrical version, the angle the 
background field makes with the grid changes as a function of 0; nonetheless, a high degree of 
symmetry is observed in the solution. We have also tested a 3D analogue of this problem on a 
200^ grid with z G [—0.5, 0.5]. Again, plots of selected variables along a horizontal line through the 
center of the blast (in the z = plane) in Figure [13] show agreement between the Cartesian and 
cylindrical versions of Athena. 

The strong S-field, /3p = 0.02 case uses the domain (x,y) G [—0.5,0.5] x [—0.5,0.5] for the 
Cartesian version, and for the cylindrical version, uses the domain (i?, 0) G [1, 2] x [2/3, 2/3], giving 
a roughly similar domain size in each case. The initial conditions consist of a circular region of 
hot gas with radius Rq = 0.125 and pressure P = 100 in an ambient medium of uniform pressure 
Pq = 1 and density po — 1. There is uniform background magnetic field of strength Bq = 10 oriented 
parallel to the x-axis. We use a computational grid of 200^ cells, third-order reconstruction and 
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the HLLD fluxes with upwind-only integration for the L/R states. The density, pressure, specific 
kinetic energy, and magnetic energy along a horizontal line through the center of the blast at 
time t = 0.02 are shown in Figure [HI with comparison to the Cartesian version of Athena. We 
have also conducted a 3D analogue of this test on a 200^ grid with z E [-0.5,0.5]. Figure [I5] 
shows a comparison of contour plots with the Cartesian version of Athena, and Figure [T6l shows a 
comparison along a horizontal line through the center of the blast (in the z = plane). 



10.7. Weber-Davis Wind 

In this problem, we invest ig!;ate a cylindrical version of the Weber-Davis wind solution as 
described in lSakurail (1 19851 ). We assume a steady, axisymmetric, 2D MHD flow with planar magnetic 
field, and a gravitational potential = —GM/R. The constants of motion are K — Pp~^, 
M RpVR, f RBr, (3 Br/{pvr), as weh as 



J 
B 



(113a) 

(113b) 
(113c) 



Here, u = {vr, — Rfl, 0) is the velocity in a frame rotating at angular velocity Q, and in this 
frame B — /3pu, so that the fluid feels no force from the magnetic field. Note that the Bernoulli 
parameter in the rotating frame includes a centrifugal potential contribution. The Alfven Mach 
number in the rotating frame is given by Ai^ = u/ca_ — l/-\/ Let Ra and pA denote the radius 
and density, respectively, at the Alfven Mach point, i.e. where J^a — 1- Then /? — 1/y^pX and 
J = R^Q,. 

Letting x = R/Ra and y = p/pA — P0^ — M~a denote the scaled radius and density, 
respectively, the Bernoulli parameter is 



GM 

~rI 




where we have defined the scaled parameters 



+ 



7-1 1 

r - - 

X 



(114) 



V = 



UJ 



Rap\GM' 

iKpY^Ra 
GM ' 
R?/^Q? 
GM ■ 



(115a) 

(115b) 
(115c) 
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Letting B = B/{GM/Ra), we have 

^dB 7j ^ ^ / {x'^-l/x^) A ^ 1 ^^^^^^ 

dx x^y^ \ (1 — y)2 / x' 

To find wind solutions, we solve the Bernoulli equation (|114l) under the constraint that the solution 
be locally flat at the slow- and fast-magnetosonic points, i.e. dB/dx — dB/dy = at (x^, yg) and 
(xj, yf). We will specify 6 and cj, and let 77, 0, Xg, y^, xj, and yj vary. Note that this becomes 
a system of si x equations in six unknowns, so if a solution exists, it must be unique. Following 
Sakurail ( 1985 ), we use the parameters 7 = 1.2, = 1. 5 and = 0.3 as for the spherically symmetric 



solar wind model of Weber and Davis. The numerical solution is given in Table [H 

We solve the problem on the domain x G [0.4, 1.8], so that the wind solution passes through 
all three critical points. We choose units such that GM — 1 and fix the initial solution at the 
inner- and outer-boundaries, evolving it long enough for equilibrium to develop. Figure [T71 shows 
the convergence of the RMS error in the Li-norm of the solution dX t — 5.0 compared to the initial 
solution. These data were computed using the Roe fluxes, second-order reconstruction, and the 
ID integrator; the results were similar for all combinations of Roe/HLLD fluxes and second-/third- 
order reconstruction. Because the 2D and 3D algorithms differ significantly from the ID version, 
especially in their treatment of magnetic fields, we present the results of the same test using these 
integrators on grids which are essentially one-dimensional, but contain a few grid cells in each 
transverse direction considered. 

Evidently, the algorithm yields second-order convergence for smooth MHD flows, and is able 
to maintain a steady magnetized solution, also conserving angular momentum as it is exchanged 
between the fluid and magnetic field. 

11. Conclusion 

We have described an adaptation of the Athena astrophysical MHD code for cyhndrical coor- 
dinates. The original Cart esian co de uses a combination of higher-order Godunov methods (based 



on the CTU algorithm of IColellal) to evolve the mass-density, momenta, and total energy, and 



constrained transport (lEvans fc Hawleyl Il988l ) to evolve the magnetic fields. We have described 
modifications to the second- and third-order reconstruction schemes, the finite- volume and finite- 
area formulations of the MHD equations, and the inclusion of geometric source terms. 

Our approach is advantageous in that it does not require modification to the majority of the 
existing code, in particular to the Riemann solvers and eigensystems. Furthermore, our approach 
to implementing cylindrical coordinates could be applied in a straightforward manner to enable 
other curvilinear coordinate systems, such as spherical coordinates, in the Athena code as well as 
other higher-order Godunov codes. 



Finally, our code and test suite are publicly available for download on the Web. The code is 
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currently being used for a variety of applications, including studies of global accretion disks, and 
we hope it will be of use to many others studying problems in astrophysical fluid dynamics. 

The authors would like to thank Jim Stone for helpful discussions about the Athena algorithms, 
and Peter Teuben for his assistance with the code package. This work was supported by grant AST 
0507315 from the National Science Foundation. 
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Fig. 8. — Convergence of the RMS error in the Li-norm for various levels of discretization of the 
rotating hydrodynamic wind test in ID, 2D and 3D. For reference, we have plotted a line of slope 
—2 (dashed) to show that the convergence is second-order in 
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Fig. 9. — For the 2D field loop advection test, we show the magnetic energy density at times 

t = and t = 2 in panels (a) and (b), respectively. Panels (c) and (d) contain magnetic field lines 
at t = and t = 2, respectively. 
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Fig. 10. — Time evolution showing dissipation of the volume- averaged magnetic energy density 
for the 2D field loop advection test. The sohd hne represents a power law fit (for t <C r) to 
the data points with residual 0.0580. 
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Fig. 11. — Contours of selected variables of the evolved state (at time t = 0.2) for the 2D MHD 
blast wave test with Bq = 1 and /3p = 0.2 using 200 x 300 grid cells, third-order reconstruction, 
HLLD fluxes, and using the cylindrical (left column) or Cartesian (right column) versions of Athena. 
Thirty equally spaced contours between the minimum and maximum are drawn in each plot. 
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Fig. 12. — Plots of selected variables along a horizontal line through the center of the blast at time 
t = 0.2 for the 2D MHD blast wave test with Bq = 1 and /3p = 0.2 using the cylindrical (circles) or 
Cartesian (solid line) versions of Athena. 
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Fig. 13. — Plots of selected variables along a horizontal line through the center of the blast at time 
t = 0.2 for the 3D MHD blast wave test with = 1 and /3p = 0.2 using the cylindrical (circles) or 
Cartesian (solid line) versions of Athena. 
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Fig. 14. — Plots of selected variables along a horizontal line through the center of the blast at time 
t = 0.02 for the 2D MHD blast wave test with = 10 and /3p = 0.02 using the cyhndrical (circles) 
or Cartesian (solid line) versions of Athena. 
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Fig. 15. — Contours of selected variables at time t = 0.02 for the 3D MHD blast wave test with 
Bo = 10 and /3p = 0.02 using 200^ grid cells and the cylindrical (top row) or Cartesian (bottom 
row) versions of Athena. Thirty equally spaced contours between the minimum and maximum are 
drawn in each plot. 
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Fig. 16. — Plots of selected variables along a horizontal line through the center of the blast at time 
t = 0.02 for the 3D MHD blast wave test with = 10 and /3p = 0.02 using the cylindrical (circles) 
or Cartesian (solid line) versions of Athena. 
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Fig. 17. — Convergence of the RMS error in the Li-norm for various levels of discretization of the 
Weber-Davis MHD wind test in ID, 2D and 3D. For reference, we have plotted a line of slope —2 
(dashed) to show that the convergence is second-order in 
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Table 1. Weber-Davis Wind Parameters 



Parameter 


Value 


7 


1.2 


e 


1.5 




0.3 


V 


2.3609 


B 


7.8745 


Xs 


0.5243 




2.4986 




1.6383 


Vf 


0.5374 



